#!/usr/bin/python 
# -*- coding: utf-8 -*- 
 
import os 
import logging 
 
reference = '/mnt/home/gongjie/workdir/c.krusei/bwa_gatk/GCF_003054445.1_ASM305444v1_genomic.fna' 
 
source_path = '/mnt/home/gongjie/workdir/c.krusei/sra/seq/'     # 此目录已经存在

names_file = source_path+'text.txt'     # 包含待处理的文件名

save_path = '/mnt/home/gongjie/workdir/c.krusei/sra/seq/'       # 此目录已经存在

# /mnt/home/gongjie/workdir/c.krusei/sra/seq/SRR6977166/

if __name__ == '__main__': 
    with open(names_file, 'r', encoding='utf-8')as f: 
        names = f.readlines() 
        
    logging.error('start....') 
    for name in names:
        name = name.strip()      # name=SRR6977166_1.clean.fastq 
        if (not name) or ('#' in name):     # 如果是空行或者使用了‘#’则该行算作不用处理的
            continue
        id = name.split('_')[0]  # id=SRR6977166    #  id   
        save_to = save_path+id+'/'
        if not os.path.exists(save_to):
            os.mkdir(save_to)
        os.system(r"bwa mem -t 4 -R '@RG\tID:HN2\tPL:Illumina\tLB:library\tSM:" + id + "' " + reference + " " + source_path + name + " > " + save_to + id + ".sam") 
        # bwa mem -t 4 -R '@RG ID:HN2 PL:Illumina LB:library SM:" + id + "' " + reference + " " + source_path + name + " > " + save_to + id + ".sam
        logging.error('111111111')
        # bwa mem -t 4 -R "@RG id:HN2 PL:Illumina LB:library SM:SRR6977166 /mnt/home/gongjie/workdir/c.krusei/bwa_gatk/GCF_003054445.1_ASM305444v1_genomic.fna /mnt/home/gongjie/workdir/c.krusei/sra/seq/SRR6977166_1.clean.fastq >SRR6977166.sam 
        os.system('samtools view -b -S ' + save_to + id + '.sam -o  ' + save_to + id + '.bam') 
        logging.error('222222222')
        os.system('samtools sort -@ 4 -m 4G -O bam -o ' + save_to + id +'.sorted.bam  ' + save_to + id + '.bam') # ????
        logging.error('333333333')
        os.system('gatk MarkDuplicates -I  ' + save_to + id + '.sorted.bam -O ' + save_to + id +'.sorted.markdup.bam -M ' + save_to + id +'.sorted.markdup_metrics.txt') 
        logging.error('444444444')
        os.system('samtools index  '+save_to + id+'.sorted.markdup.bam') 
        logging.error('5555555555')
        os.system('srun gatk HaplotypeCaller -R  '+reference+' -I  '+save_to + id+'.sorted.markdup.bam -O  '+save_to + id+'.g.vcf -ERC GVCF') 
        logging.error('66666666666')
        os.system('gatk GenotypeGVCFs -R  '+reference+' -V  '+save_to + id+'.g.vcf -O  '+save_to + id+'.HC.vcf') 
        logging.error('77777777777')
        os.system('gatk SelectVariants -R  '+reference+' -select-type SNP -V  '+save_to + id+'.HC.vcf -O  '+save_to + id+'.HC.SNP.vcf') 
        logging.error('88888888888')
        os.system('gatk VariantFiltration -R  '+reference+' -V  '+save_to + id+'.HC.SNP.vcf --filter "QD <2.0 || MQ <40.0 || FS >60.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0 || HaplotypeScore >13.0" --filter-name "LowQual" -O  '+save_to + id+'.HC.SNPs.filter.vcf') 
        logging.error('99999999999')
        os.system('grep -v "LowQual"  '+save_to + id+'.HC.SNPs.filter.vcf > '+save_to + id+'.HC.SNPs.pass.vcf') 
        #save_to + id+'.HC.SNPs.filter.vcf
    logging.error('end....') 


 
 